#include <fstream>
#include <iostream>
#include <iomanip>
#include "TH1.h"
#include "TF1.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMath.h"
#include "TString.h"
#include "TFeldmanCousins.h"

using namespace std;

TFile *ff;

double AA = 40.; //Ar=40 Xe=131.3
double Day = 100.; //day
double Volume =100.; //kg
double Nobserved   = 0.; //N of Observed Event
double Nbackground = 0.; //N of Expected BG Event

double Ew = 19.5e-3; // (keV)
int    ThrePE = 10; //Thereshold of detected photon
double Efficiency = 0.5; // efficiency

//------------------------------------------- Parameters
// A : atmic mass
// Mt : target mass
// Mwimp : wimp mass
// sigmawn : Xsec for wimp and nucleon
double c = 3.0e8; // [m/s]
double ck = 3.0e5; // [km/s]
double rho0 = 0.3; //[GeV/c2/cm3] WIMP energy density
double Mwimp; // [GeV/c2] WIMP mass
double Mn = 0.9315; // [GeV/c2] average nucleon mass
double v0 = 220;    // [km/s] Sun's velocity from Galactic center
double vesc = 544;  // [km/s] Galactic escape velocity
//double vesc = 600;  // [km/s] Galactic escape velocity
double vear = 232;  // [km/s] Average earth velocity from Galactic center
double NA = 6.0221415e23;  // Avogadro's number
double pi = TMath::Pi();

double c1 = 0.751;      //factor change for each months
double c2 = 0.561;      //factor change for each months
//-------------------------------------------------------------

//double sigma0(double A, double Mt, double Mw, double sigmawn); 
double Mr(double M1, double M2); // reduced mass
double F(double Q, double A, double mT);

//-----
void SetWIMPMass(double WM){
  Mwimp = WM; //[GeV/c2] WIMP mass
};

//-----

double dRdE(double ER ,double A , double sigmawn){
  // argon 
  double Mt = A * 0.932; // GeV/c2
  double sigmaA = pow(A,2) * pow( Mr(Mwimp,Mt)/Mr(Mwimp,Mn),2 ) * sigmawn;
  double vmin = sqrt( Mt*ER / (2*Mr(Mwimp,Mt)) );
  double Q = sqrt(2*Mt*ER)/c; // keV s/m

  double k0 = pow(pi*pow(v0,2),3/2);
  double kf = k0 * ( TMath::Erf(vesc/v0) - 2/sqrt(pi)*vesc/v0*exp( -pow(vesc/v0,2) ) );
  //  double R0 = 2 *NA * rho0/sqrt(pi)/A/Mwimp * sigmaA * v0;
  double R0 = (361/Mwimp/Mt) * (sigmaA) * (rho0/0.3) * (v0/220);   // GeV , pb , GeV/cm3 , km/s  ## Nishimura ph.D ##

  double E0 = Mwimp*pow((v0/ck),2)/2.;         // GeV*(km/s)^2   E(100GeV)=mc^2
  double rm = 4*Mwimp*Mt/pow(Mwimp+Mt,2); 

  //  double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf((vmin+vear)/v0) - TMath::Erf((vmin-vear)/v0) );
  double Term1 = c1 * exp(-c2 * ER / E0 / rm);  // year average
  double Term2 = exp( -pow(vesc/v0,2) );
 
  //***1
  //  double drate =  R0/E0/rm * exp(-ER/E0/rm) *1e-6;  // vesc = inf , vear = 0
  //***2
  double drate =  R0/E0/rm * Term1 *1e-6;  // vesc = inf , vear = vear
  //***3
  //double drate = k0/kf * R0/E0/rm * ( Term1 - Term2 ) *1e-6; //  vesc = vesc , vear = vear 

    //  return  R0/E0/rm * exp(-ER/E0/rm) *1e-6;
  return drate;//100kg*100day
}


double dRdE2(double *x ,double *par){
  // argon 
  double ER = x[0]*1e-6; // keV
  double A = par[0];  // atomic mass
  double Mt = A * 0.932; // GeV/c2
  double sigmawn = par[1]; // pb

  double sigmaA = pow(A,2) * pow( Mr(Mwimp,Mt)/Mr(Mwimp,Mn),2 ) * sigmawn;
  double vmin = sqrt( Mt*ER / (2*Mr(Mwimp,Mt)) );
  double Q = sqrt(2*Mt*ER)/c; // keV s/m

  double k0 = pow(pi*pow(v0,2),3/2);
  double kf = k0 * ( TMath::Erf(vesc/v0) - 2/sqrt(pi)*vesc/v0*exp( -pow(vesc/v0,2) ) );
  //  double R0 = 2 *NA * rho0/sqrt(pi)/A/Mwimp * sigmaA * v0;
  double R0 = (361/Mwimp/Mt) * (sigmaA) * (rho0/0.3) * (v0/220);   // GeV , pb , GeV/cm3 , km/s

  double E0 = Mwimp*pow((v0/ck),2)/2.;         // GeV*(km/s)^2   E(100GeV)=mc^2
  double rm = 4*Mwimp*Mt/pow(Mwimp+Mt,2); 

  //  double Term1 = sqrt(pi)*v0/(4*vear) * ( TMath::Erf((vmin+vear)/v0) - TMath::Erf((vmin-vear)/v0) );
  double Term1 = c1 * exp(-c2 * ER / E0 / rm);
  double Term2 = exp( -pow(vesc/v0,2) );
  
  //  double dRdE = k0/kf * R0/E0/rm * ( Term1 - Term2 ); //* pow(F(Q,A,Mt),2);
  double drate =  R0/E0/rm * exp(-ER/E0/rm) *1e-6;

  return drate; //100kg*100day
}

// reduce mass
double Mr(double M1, double M2){
  return M1*M2/(M1+M2);
}

//form factor  Helm form factor (CDMS thesis)           
// ER -> recoil E , mT -> target mass (keV)
double F(double ER, double A, double mT){
  int option=0;
  
  double Q = sqrt(2*mT*ER)/c; // keV s/m                                              

  double aT = 0.52                    *pow(10,-15); // m                              
  double cT = (1.23*pow(A,1/3.) - 0.6)*pow(10,-15); // m                              
  double sT = 0.9                     *pow(10,-15); // m                              

  double rT = sqrt( pow(cT,2)+7/3.*pow(pi*aT,2)-5*pow(sT,2) ); // [m] effective nuclear radius  

  double hbar = 6.582e-19; // keVs                                                    
  double xq = Q*rT/hbar; // dimentionless                                             
  
  double formfac = 3*(sin(xq) - xq*cos(xq))/pow(xq,3) * exp( -pow(Q*sT/hbar,2)/2. );       
  if(option==1) formfac = 3*sqrt(pi/(2*xq))*TMath::BesselJ1(xq)/xq * exp( -pow(Q*sT/hbar,2)/2. );  // use TMath BesselJ1 if you want 

 return formfac;
}

//-----------------                                                                   
// test form factor                                                                   
//-----------------                                                                   
double FormFac(double *x, double *par){
  double ER1 = x[0];
  double A1 = par[0];
  double mT1 = par[1];
  return pow(F(ER1,A1,mT1),2);
}

// for plot  -> dRdE * FormFactor

double dForm(double *x, double *par){
  double er = x[0];
  double A  = par[0];
  double sigmawn  = par[1];
  double mT = par[2];
  
  double er2=er * 1e-6;
  double drF = Day*Volume*dRdE(er2,A,sigmawn) * pow(F(er,A,mT),2); //Changed By Sugita For day & volume

  //  double drF = dRdE(er2,A,sigmawn);
  //return drF;
  
  if(drF<0){return 0;}
  else     {return drF;}
}

Double_t myfunc(double *x,double *par){

  Double_t fff1=(TMath::Erf((x[0]-par[0])/par[1])+1.0)*0.5*par[2];
  Double_t fff2=(TMath::Erf((x[0]-par[3])/par[4])+1.0)*0.5*(1.0-par[2]);

  if( (fff1+fff2)<0.0001 ) return 0;
  //if( (fff1+fff2)<0.005 ) return 0;
  else                    return fff1+fff2;
}

double dForm_eff(double *x, double *par){
  //return  dForm(x,par)*myfunc(x,&par[3]);
  return  dForm(x,par)*myfunc(x,&par[3])*0.5;
}

void DMeff(double WM,int IsSave=0){

  SetWIMPMass(WM);

  int RR=20;
  int CONFIG=2;
  float REF=0.9;

  ff=new TFile(Form("./rootfile/R%d_config%d_Alpha0.18_ref%.1f.root",RR,CONFIG,REF));
  TH1F* heff=(TH1F*)ff->Get("hEffvsEr");

  TF1 *func=new TF1("func",myfunc,0,1000,5);
  func->SetParLimits(0,0,100.);
  func->SetParLimits(1,0.,100.);
  func->SetParLimits(2,0.1,1.);
  func->SetParLimits(3,10.,100);
  func->SetParLimits(4,10.,100.);
  heff->Fit("func","Q");

  double ParSigma = 1e-44; //Asume sigma 
  TF1 *f0 = new TF1("f0", dForm, 0, 1000, 3);
  TF1 *f1 = new TF1("f1", dForm_eff, 0, 1000, 9); 
  f0->SetParameters(AA,ParSigma*1e36,AA*0.932e6);   // A->Ar(40), sigma , mT =A*0.932(keV)
  f1->SetParameters(AA,ParSigma*1e36,AA*0.932e6,
		    func->GetParameter(0),
		    func->GetParameter(1),
		    func->GetParameter(2),
		    func->GetParameter(3),
		    func->GetParameter(4),
		    func->GetParameter(5)
		    ); 

  double integral     = f0->Integral(0,1000);
  double integral_eff = f1->Integral(0,1000);

  double ratio = integral_eff / integral;

  cout<< integral_eff<<" / "<<integral<<" = "<<ratio<<endl;

  TCanvas *can = new TCanvas("can", "can");
  //  can->SetLogx();
  can->SetLogy();
  can->SetGrid(1);
  TH1 *fr=gPad->DrawFrame(10,1e-5,200,10);
  fr->SetTitle(Form("Wimp Mass = %.0f GeV, #sigma = %.1e cm^{2}",WM,ParSigma));
  fr->SetXTitle("Recoil Energy (keV)");
  fr->SetYTitle(Form("Count /(keV*%.0fkg*%.0fday)",Volume,Day));
  f0->SetLineColor(1);
  f1->SetLineColor(4);
  fr->Draw();
  f0->Draw("same");
  f1->Draw("same");

  if(IsSave){  
    char name[128];
    sprintf( name, "./txt/R%d_config%d_Alpha0.18_ref%.2f.txt",RR,CONFIG,REF );
    ofstream fout( name ,std::ios::out | std::ios::app);
    fout<<WM<<" "<<ratio<<endl;
  }
  cout<<"WIMP = "<<WM<<" GeV, Acceptance = "<<ratio*100<<"%"<<endl;

}
